# t-test with complete output
# df for two sample t test with unequal variances
# function written by Guadalupe Tuñón
###################################################

ttest <- function(y, x, two.tailed=TRUE){ 
  
  # Calculating difference in means
  mean1 <- mean(y[x==1], na.rm=T)
  mean0 <- mean(y[x==0], na.rm=T)
  diff <- mean1 - mean0
  
  # Calculating SE of the difference
  N1 <- length(na.omit(y[x==1]))
  N0 <- length(na.omit(y[x==0]))
  var1 <- var(y[x==1],na.rm=T)
  var0 <- var(y[x==0],na.rm=T)
  varN1 <- var1/N1
  varN0 <- var0/N0
  se.diff <- sqrt(varN1 + varN0)
  
  # T-statistic
  t <- diff/se.diff
  
  # Degrees of freedom
  df.num <- ((varN1 + varN0)^2)
  df.den <- (varN1^2)/(N1-1) + (varN0^2)/(N0-1)
  df <- df.num/df.den
  
  # P-value
  if(two.tailed==TRUE){
    p <- 2 * pt(abs(t), df, lower.tail=F)
  }
  
  if(two.tailed==FALSE){
    p <- pt(t, df, lower.tail=F)
  }
  
  N <- N1 + N0
  
  # Preparing output
  res <- c(mean1, mean0, diff, se.diff, 
           t, N, df, p)
  names(res) <- c("Mean 1", "Mean 0", "Difference", 
                  "SE_Diff","t-stat", "N", "df", "p-value")
  
  return(c(res))
}



anderson.rubin.ci=function(ivmodel,conflevel=.95){
  y=ivmodel$y
  n=length(y);
  if(is.null(ivmodel$x)==TRUE){
    print("Refit ivmodel with x=TRUE option in the ivreg function.")
    stop()
  }
  regressors=ivmodel$x$regressors
  instruments=ivmodel$x$instruments
  # Use notation in Davidson and MacKinnon, 2011
  # Figure out based on ivmodel fit from ivreg, the elements in Davidson and MacKinnon 
  W=instruments;
  # Figure out which columns in regressors and instruments are the same
  regressors.same.vec=rep(0,ncol(regressors))
  instruments.columns.same=rep(0,ncol(regressors)-1)
  count=1
  for(i in 1:ncol(regressors)){
    tempmat=regressors[,i]-instruments
    tempsum=apply(abs(tempmat),2,sum)
    regressors.same.vec[i]=sum(tempsum==0)>0
    if(sum(tempsum==0)>0){
      instruments.columns.same[count]=i
      count=count+1
    }
  }
  
  y2column=which(regressors.same.vec==0)
  y2=as.matrix(regressors[,y2column])
  Z=as.matrix(regressors[,-y2column])
  W2=as.matrix(instruments[,-instruments.columns.same])
  
  l=ncol(W);
  k=ncol(Z);
  q=qf(conflevel,l-k,n-l)
  cval=q*(l-k)/(n-l)
  y1=matrix(y,ncol=1)
  
  y2hat.given.W=fitted(lm(y2~W-1))
  y2hat.given.Z=fitted(lm(y2~Z))
  y1hat.given.W=fitted(lm(y1~W-1))
  y1hat.given.Z=fitted(lm(y1~Z))
  coef.beta0sq=cval*sum(y2^2)-(cval+1)*sum(y2*y2hat.given.W)+sum(y2*y2hat.given.Z)
  coef.beta0=-2*cval*sum(y1*y2)+2*(cval+1)*sum(y1*y2hat.given.W)-2*sum(y1*y2hat.given.Z)
  coef.constant=cval*sum(y1^2)-(cval+1)*sum(y1*y1hat.given.W)+sum(y1*y1hat.given.Z)
  D=coef.beta0^2-4*coef.constant*coef.beta0sq
  if(coef.beta0sq==0){
    if(coef.beta0>0){
      ci=c("[",-coef.constant/coef.beta0,",Infinity)")
    }
    if(coef.beta0<0){
      ci=c("(-Infinity,",-coef.constant/coef.beta0,"]");
    }
    if(coef.beta0==0){
      if(coef.constant>=0){
        ci="Whole Real Line"
      }
      if(coef.constant<0){
        ci="Empty Set"
      }
    }
  }
  if(coef.beta0sq!=0){
    if(D<0){
      if(coef.beta0sq>0){
        ci="Whole Real Line"
      }
      if(coef.beta0sq<0){
        ci="Empty Set"
      }
    }
    if(D>0){
      # Roots of quadratic equation
      # Roots of quadratic equation
      root1=(-coef.beta0+sqrt(D))/(2*coef.beta0sq)
      root2=(-coef.beta0-sqrt(D))/(2*coef.beta0sq)
      upper.root=max(root1,root2)
      lower.root=min(root1,root2)
      if(coef.beta0sq<0){
        ci=paste("[",lower.root,",",upper.root,"]")
      }
      if(coef.beta0sq>0){
        ci= paste("(-Infinity,",lower.root,"] union [",upper.root,",Infinity)")
      }
    }
    if(D==0){
      ci="Whole Real Line"
    }
  }
  list(confidence.interval=ci)
}

ARsensitivity.ci=function(ivmodel, Delta=NULL, conflevel=.95){
  if(is.null(ivmodel$x)==TRUE){
    print("Refit ivmodel with x=TRUE option in the ivreg function.")
    stop()
  }
  regressors=ivmodel$x$regressors
  instruments=ivmodel$x$instruments
  
  # Figure out which columns in regressors and instruments are the same
  regressors.same=rep(0,ncol(regressors))
  instruments.same=rep(0,ncol(instruments))
  
  for(i in 1:ncol(regressors)){
    tempmat=regressors[,i]-instruments
    tempsum=apply(abs(tempmat),2,sum)
    regressors.same[i]=sum(tempsum==0)>0
    if(sum(tempsum==0)>0){
      instruments.same[which(tempsum==0)]=i
    }
  }
  
  if(sum(regressors.same==0)!=1){
    print("Input exact one exposure for the sensitivity method.")
    stop()
  }
  if(sum(instruments.same==0)!=1){
    print("Input exact one IV for the sensitivity method.")
    stop()
  }
  
  D=as.matrix(regressors[, which(regressors.same==0)])
  X=as.matrix(regressors[, -which(regressors.same==0)])
  Z=as.matrix(instruments[, which(instruments.same==0)])
  Y=matrix(ivmodel$y, ncol=1)
  W=instruments
  Zhat.given.Xperp=residuals(lm(Z~X-1))
  
  if(!is.null(Delta)){
    if(is.numeric(Delta) & length(Delta)==2){
      ncp=max(Delta^2)*sum(Zhat.given.Xperp^2)
    }else{
      print("Wrong input of the sensitivity range.")
      stop()
    }
  }else{
    ncp=0
  }
  
  n=nrow(Y)
  k=ncol(X);
  q=qf(conflevel, df1=1, df2=n-k-1, ncp=ncp)
  cval=q/(n-k-1)
  
  Yhat.given.W=fitted(lm(Y~W-1))
  Yhat.given.X=fitted(lm(Y~X-1))
  Dhat.given.W=fitted(lm(D~W-1))
  Dhat.given.X=fitted(lm(D~X-1))
  
  coef.beta0sq=cval*sum(D^2)-(cval+1)*sum(D*Dhat.given.W)+sum(D*Dhat.given.X)
  coef.beta0=-2*cval*sum(D*Y)+2*(cval+1)*sum(Y*Dhat.given.W)-2*sum(Y*Dhat.given.X)
  coef.constant=cval*sum(Y^2)-(cval+1)*sum(Y*Yhat.given.W)+sum(Y*Yhat.given.X)
  D=coef.beta0^2-4*coef.constant*coef.beta0sq
  
  ci=matrix(NA, ncol=2)
  colnames(ci)<-c("lower", "upper")
  
  if(coef.beta0sq==0){
    if(coef.beta0>0){
      info=c("[",-coef.constant/coef.beta0,",Infinity)")
      ci[1,]=c(-coef.constant/coef.beta0, Inf)
      type=1
    }
    if(coef.beta0<0){
      info=c("(-Infinity,",-coef.constant/coef.beta0,"]");
      ci[1,]=c(-Inf, -coef.constant/coef.beta0)
      type=1
    }
    if(coef.beta0==0){
      if(coef.constant>=0){
        info="Whole Real Line"
        ci[1,]=c(-Inf, Inf)
        type=1
      }
      if(coef.constant<0){
        info="Empty Set"
        type=2
      }
    }
  }
  if(coef.beta0sq!=0){
    if(D<0){
      if(coef.beta0sq>0){
        info="Whole Real Line"
        ci[1,]=c(-Inf, Inf)
        type=1
      }
      if(coef.beta0sq<0){
        info="Empty Set"
        type=2
      }
    }
    if(D>0){
      # Roots of quadratic equation
      # Roots of quadratic equation
      root1=(-coef.beta0+sqrt(D))/(2*coef.beta0sq)
      root2=(-coef.beta0-sqrt(D))/(2*coef.beta0sq)
      upper.root=max(root1,root2)
      lower.root=min(root1,root2)
      if(coef.beta0sq<0){
        info=paste("[",lower.root,",",upper.root,"]")
        ci[1, ]=c(lower.root, upper.root)
        type=0
      }
      if(coef.beta0sq>0){
        info= paste("(-Infinity,",lower.root,"] union [",upper.root,",Infinity)")
        ci[1, ]=c(-Inf, lower.root)
        ci<-rbind(ci, c(upper.root, Inf))
        type=1
      }
    }
    if(D==0){
      info="Whole Real Line"
      ci[1, ]=c(-Inf, Inf)
      type=1
    }
  }
  list(confidence.interval=ci, printinfo=info, ci.type=type)
}

ARsensitivity.power <-function(n,k,lambda,gamma,var.z,sigma1,sigma2,rho,alpha=.05,
                               Delta=NULL, delta=NULL){
  if(!is.null(Delta)){
    if(is.numeric(Delta) & length(Delta)==2){
      Delta <- sort(Delta)
      if(!is.null(delta)){
        if(is.numeric(delta) & length(delta)==1){
          ncp1 = (lambda*gamma+delta*sigma1)^2*n*var.z/
            (sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
          ncp2 = max(Delta^2)*n*var.z
        }else{
          print("Wrong input of the true sensitivity parameter")
          stop()
        }	
      }else{
        if((Delta[1] < -gamma*lambda/sigma1) & (Delta[2] > -gamma*lambda/sigma1)){
          ncp1=0
        }else{
          ncp1=min((lambda*gamma+Delta[1]*sigma1)^2,(lambda*gamma+Delta[2]*sigma1)^2)*
            n*var.z/(sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
        }
        ncp2 = max(Delta^2)*n*var.z
      }
    }else{
      print("Wrong input of the sensitivity range.")
      stop()
    }
  }else{
    if(!is.null(delta)){
      print("Sensitivity range missing")
    }else{
      ncp1 = lambda^2*gamma^2*n*var.z/
        (sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
      ncp2 = 0
    }
  }
  
  temp = qf(1-alpha, df1=1, df2=n-k-1, ncp=ncp2)  
  power = 1-pf(temp, df1=1, df2=n-k-1, ncp=ncp1)
  return(power)
}

ARsensitivity.size <-function(power,k,lambda,gamma,var.z,sigma1,sigma2,rho,alpha=.05,
                              Delta=NULL, delta=NULL){
  if(!is.null(Delta)){
    if(is.numeric(Delta) & length(Delta)==2){
      Delta <- sort(Delta)
      if(!is.null(delta)){
        if(is.numeric(delta) & length(delta)==1){
          ncp1 = (lambda*gamma+delta*sigma1)^2*var.z/
            (sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
          ncp2 = max(Delta^2)*var.z
        }else{
          print("Wrong input of the true sensitivity parameter")
          stop()
        }	
      }else{
        if((Delta[1] < -gamma*lambda/sigma1) & (Delta[2] > -gamma*lambda/sigma1)){
          ncp1=0
        }else{
          ncp1=min((lambda*gamma+Delta[1]*sigma1)^2,(lambda*gamma+Delta[2]*sigma1)^2)*
            var.z/(sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
        }
        ncp2 = max(Delta^2)*var.z
      }
    }else{
      print("Wrong input of the sensitivity range.")
      stop()
    }
  }else{
    if(!is.null(delta)){
      print("Sensitivity range missing")
      stop()
    }else{
      ncp1 = lambda^2*gamma^2*var.z/
        (sigma1^2+2*rho*sigma1*sigma2*lambda+sigma2^2*lambda^2)
      ncp2 = 0
    }
  }
  
  if(ncp1<=ncp2){
    print("Sensitivity range too large")
    stop()
  }else{
    oldn<-k+2
    state<-1
    while(state){
      temp = qf(1-alpha, df1=1, df2=oldn-k-1, ncp=ncp2*oldn)  
      temppower = 1-pf(temp, df1=1, df2=oldn-k-1, ncp=ncp1*oldn)  
      if(temppower < power){
        oldn <- oldn*2
      }else{
        state <- 0
      }
    }
    
    lower <- oldn%/%2
    upper <- oldn
    while((upper-lower)>2){
      new <- (upper+lower)%/%2
      temp = qf(1-alpha, df1=1, df2=oldn-k-1, ncp=ncp2*new)  
      temppower = 1-pf(temp, df1=1, df2=oldn-k-1, ncp=ncp1*new)  
      if(temppower < power){
        lower <- new
      }else{
        upper <- new
      }
    }
  }
  return(upper)
}


cluster.robust.se <-
  function(ivmodel,clusterid){
    print("Cluster Robust Standard Errors");
    clx(ivmodel,clusterid);
  }


clx <-
  function(fm,cluster){
    dfcw=1
    M <- length(unique(cluster))
    N <- length(cluster)
    dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
    u <- apply(estfun(fm),2,
               function(x) tapply(x, cluster, sum))
    vcovCL <- dfc*sandwich(fm, meat.=crossprod(u)/N)*dfcw
    vcovCL
    coeftest(fm, vcovCL) 
  }

power.iv <-
  function(n,lambda,gamma,var.z,sigmau,sigmav,rho,alpha=.05){
    Lambda=(lambda^2)/((sigmau/sigmav)^2+2*rho*(sigmau/sigmav)*lambda+lambda^2);
    ncp=gamma^2*(n*var.z)*Lambda/sigmav^2;
    fquantile=qf(1-alpha,1,n-2)
    power=1-pf(fquantile,1,n-2,ncp);
    list(power=power);
  }

robust.se <-
  function(ivmodel){
    robustcov=sandwich(ivmodel);
    print("Robust Standard Errors")
    coeftest(ivmodel,robustcov);
  }

sigmav.func <-
  function(prob.d1.given.z1,prob.d1.given.z0,prob.z1){
    probvector=c((1-prob.z1)*prob.d1.given.z0,(1-prob.z1)*(1-prob.d1.given.z0),prob.z1*prob.d1.given.z1,prob.z1*(1-prob.d1.given.z1));
    value.vector=c(1-prob.d1.given.z0,-prob.d1.given.z0,1-prob.d1.given.z1,-prob.d1.given.z1);
    expected.value=sum(probvector*value.vector);
    variance=sum(probvector*value.vector^2)-expected.value^2;
    list(sigmav=sqrt(variance))
  }





